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Simulating Bone Loss in Microgravity Using Mathematical 
Formulations of Bone Remodeling 

James A. Pennline 

National Aeronautics and Space Administration 
Glenn Research Center 
Cleveland, Ohio 44135 

Abstract 

Most mathematical models of bone remodeling are used to simulate a specific bone disease, by 
disrupting the steady state or balance in the normal remodeling process, and to simulate a therapeutic 
strategy. In this work, the ability of a mathematical model of bone remodeling to simulate bone loss as a 
function of time under the conditions of microgravity is investigated. The model is formed by combining 
a previously developed set of biochemical, cellular dynamics, and mechanical stimulus equations in the 
literature with two newly proposed equations; one governing the rate of change of the area of cortical 
bone tissue in a cross section of a cylindrical section of bone and one governing the rate of change of 
calcium in the bone fluid. The mechanical stimulus comes from a simple model of stress due to a 
compressive force on a cylindrical section of bone which can be reduced to zero to mimic the effects of 
skeletal unloading in microgravity. The complete set of equations formed is a system of first order 
ordinary differential equations. The results of selected simulations are displayed and discussed. 
Limitations and deficiencies of the model are also discussed as well as suggestions for further research. 

1.0 Introduction 

Under the conditions of microgravity, astronauts lose bone mass and calcium (Ref. 1). In fact, bone can 
atrophy in space at a rate of about 1 percent a month, and some models suggest that a 40 to 60 percent total 
loss could eventually be reached (Ref. 13). Data from Buckley (Ref. 1) indicates that the greatest loss occurs 
in the lower extremities which experience the higher skeletal loadings during normal activity on earth. 
National Aeronautics and Space Administration (NASA) has a specific interest in understanding the 
mechanisms for bone loss during space flight for several reasons. Bone recovery has proven problematic. 
Based on a recovery model developed from medical data on astronauts who had flown on long duration 
missions (4 to 6 months), it could take as much as 9 months to recover 50 percent of the bone loss in some 
skeletal sites (Ref. 28). Other health hazards instigated by skeletal bone loss are accumulations of excess 
minerals in tissue such as the kidney, and increased risk of fracture. NASA is interested in developing 
interventions to sustain and protect the well being of its astronauts. Consequently, understanding bone loss 
is one of the research efforts being conducted for development of The Digital Astronaut. The Digital 
Astronaut is a an integrated, modular modeling and database system intended to support biomedical 
research to help identify and interpret medical and physiological research, to determine the effectiveness of 
specific human countermeasures to meet health and performance goals on exploration missions, and to 
evaluate medical interventions during mission emergencies, accidents, and illnesses (Ref. 24). Of particular 
concern is the effect on bone health from a longer duration mission such as a journey to Mars or the 
development of a Moon base. Existing measures of scheduled exercise with specially designed equipment 
may not be sufficient for preventing bone loss. 

Bone remodeling is the process by which bone is removed and replaced on the same surface site by 
highly mediated and coupled action of bone cells. It’s the physiological mechanism for enabling the 
normal maintenance of bone in the adult and is achieved by a balance between the processes of bone 
formation and bone resorption. These processes can be influenced by skeletal loading, endocrine 
regulation, and local biochemical mediators. Bone loss experienced by astronauts is thought to be caused 
by a disruption in the resorption-formation balance in favor of resorption, triggered by skeletal unloading. 
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Gravity is somehow converted from a mechanical signal to a chemical signal. Although much is known 
about the chemical signals, the process by which cells sense mechanical stimuli and then translate the 
information to a signal (mechanotransduction) that causes response within the cell or another cell, is not 
fully understood (Refs. 13 and 2). 

Significant knowledge has been gained in understanding the biomechanical and molecular regulation 
of bone remodeling (Ref. 20). Work done in biomechanics has used mathematical formulations to 
describe how bone responds to stress and strain under loading. For example, in the study of fracture 
mechanics, mechanical and material properties are used to study how bones deform under loads. At the 
molecular and cellular level, building mathematical models of cellular dynamics has become an 
increasingly active research area in the last decade (Refs. 4, 5, 7, 8, 9, 14, 17, 19, 23, and 25). The cells 
that carry out the remodeling process are the osteoclasts that resorb bone by secretion of acid and 
proteases, and the osteoblasts that form bone by forming an initial collagen matrix (osteoid) and then 
mineralizing the collagen (Ref. 15). Osteocytes, matured osteoblasts encased in calcified bone matrix, are 
the cells that have been suggested to be the stimulus sensors or receptors of the stimulus signal. This 
suggestion is consistent with histologic and physiologic data (Ref. 2). One of the more descriptive 
models, which attempts to accurately detail the biochemical process, is the one by Lemaire et al. (Ref. 7). 
That model has also been used by other researchers as a base on which to build an extended model 
(Refs. 9, 14, and 17). 

The majority of models in the literature lack some elements that prevent them from being used to model 
and simulate the effects of microgravity. In particular, there is a lack of a sufficient way to quantify bone 
loss or gain. The model in (Ref. 4) describes the change in bone density, as a percent, during a bone 
remodeling units cycle as a function of the change in the number of cells. Models in (Refs. 14 and 17) use 
similar representations of changes in relative mass. None of the others specify an explicit measure of bone 
loss or gain. Models also need to link the effect of the dynamics of loads on bone to the metabolic process. 
The model in (Ref. 9) is the first one this author has seen that attempts to combine how bone adapts to both 
cellular interactions and mechanical stimulus. It formulates the change in the intracortical radius of a 
cylindrical section of cortical bone under a force-induced stimulus that influences the cellular dynamic 
effects. 

Most of the current cellular models focus on simulating the cause of specific bone diseases that occur 
on earth and the result of possible therapeutic strategies. The objective of this work is to investigate the 
applicability of some of the cellular dynamic models to NASA’s interest in developing a model to help 
simulate the effects of microgravity and to provide quantitative analysis of bone mass and calcium levels 
to Digital Astronaut (DA). To do so we describe a model that uses the cellular dynamics from (Ref. 7), 
the force induced stimulus function from (Ref. 9), a new equation quantifying bone loss or gain, and an 
additional mass balance equation for calcium in bone fluid. In particular, the equation in (Ref. 9) that 
governs the changes in the intracortical thickness of a cylindrical section of cortical bone is replaced with 
an equation that governs the changes in the area of hard bone tissue in a cross section. The new equation 
can be related to the rate of change in the complement of bone porosity as defined in (Refs. 1 1 and 12). 
Section 2 gives a brief description of bone anatomy and the remodeling process. In Section 3, the physical 
part of the model is described and in Section 4, the model equations and variables are defined. In Section 
5, the results of example simulations are discussed. Conclusions and recommendations for future research 
are made in Sections 6. 


2.0 Bone Remodeling Process 

The skeletal makeup includes cortical bone, the compact bone that forms the hard outside shell that 
encloses bone marrow, and trabecular (cancellous) bone, the spongy interior tissue in bone marrow. Trabecular 
bone is made up of a network of rod and plate like elements. Cortical bone is organized into Haversian 
systems, consisting of concentric layers of bone and a central canal (Haversian canal) which supplies blood. 
The boundaries between Haversian systems are referred to as cement lines. Within bone, osteocytes, cells 
derived from the bone forming cells, form what is understood to be a signaling network (Ref. 31). 
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Figure 1. — Compact bone and spongy (cancellous bone). 


Bone remodeling, the physiological mechanism for maintenance, renewal, and repair in the adult 
skeleton, is done threw the replacement of bone in units by the coupled action of bone cells on the same 
bone surface. The bone resorbing cells, osteoclasts, remove old or damaged bone. The bone forming cells, 
osteoblasts, form an initial collagen matrix and then mineralize the collagen. The replacement unit, often 
referred to as a basic multicellular unit (BMU) or bone remodeling unit (BRU), differs between trabecular 
bone and cortical bone. In trabecular bone the structural unit is a packet ideally shaped like a shallow 
crescent (hemi-osteon) on the surface of a rod or plate like element. In cortical bone the structural unit is a 
single Haversian system (Osteon) shaped like a cylinder and is referred to as a cutting cone while forming 
(Ref. 32). 

The remodeling cycle of a BRU consists of 5 phases, activation, resorption, reversal, formation, and 
quiescence. Activation, the conversion of a small area of surface from quiescence to activity with respect 
to remodeling, requires recruitment of osteoclasts, a means for them to gain access to bone, and a 
mechanism for their attraction and attachment to the surface. Because of variations in age, sex, race, and 
metabolic state, and systematic differences between different bones and between different surfaces in the 
same bone, activation occurs sometimes at random and sometimes in response to focal events in the bone 
structure or biomechanical stimulus. In the normal adult skeleton it occurs about once every 10 sec 
(Ref. 15). Following activation, osteoclasts, proliferated from precursor cells, begin to erode a pit referred 
to as a Howship’s lacuna in trabecular bone or burrow a tunnel called a cutting cone in cortical bone. (The 
pit characteristics were described in the previous paragraph.) The completion of resorption can take from 
one to three weeks. Reversal, the time interval between completion of resorption and beginning of 
formation at the same location, can last from one to two weeks. Formation by osteoblasts, derived from a 
different precursor cell, involves collagen matrix formation and mineralization. It lasts about three 
months, but it can take an additional several months for newly formed bone to become mature. Thus, the 
complete cycle is on the average about four months (Ref. 15). There are an estimated 10 5 to 10 6 active 
BRUs in the entire body at any one time (Ref. 33). 

The processes of remodeling are controlled or influenced by hormones, protein receptors, ligands, 
other biochemicals, and mechanical stimulus. A complete description of the entire process of endocrine 
regulation, biochemical mediation, and stimulation by skeletal loading involved in bone remodeling is 
beyond the scope of this report. We shall briefly discuss the main substances involved in the cellular 
process that appear in the model equations of Section 4. The reader can refer to a list of their definitions 
in Appendix B. Early in the remodeling cycle, surface bound molecules, referred to as Receptor Activator 
for Nuclear Factor kB ligand (RANKL), are expressed on the surface of osteoblasts while the membrane 
protein Receptor Activator for Nuclear Factor kB (RANK) is expressed on the surface of preosteoclasts. 
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The binding of the receptor RANK with the ligand RANKL causes the derivation of active osteoclasts 
from preosteoclasts (Ref. 38). Parathyroid hormone (PTH), secreted by the parathyroid gland, is involved 
in promoting bone resorption by binding to osteoblasts and stimulating them to increase expression of 
RANKL. Thus PTH indirectly stimulates an increase in osteoclasts. PTH can have opposite effects on 
bone mass depending on whether it is released into the plasma continuously or on a pulsatory schedule 
(Ref. 5). Osteoprotegerin (OPG) is released by osteoblasts and acts as a decoy receptor for RANKL which 
blocks RANK-RANKL binding and inhibits preosteoclasts from maturing into active osteoclasts 
(Ref. 38). Transforming Growth Factor Beta (TGF-P) is released by osteoclasts during bone resorption. It 
has the ability to stimulate an increase of and activity in the osteoblast population, but it can also inhibit 
final maturation into active osteoblasts (Refs. 35 and 36). It’s also thought to induce osteoclast apoptosis 
(Ref. 37). Nitric oxide, a gaseous signaling molecule, produced by osteoblasts directly regulates 
osteoclastic activity (Ref. 34). Prostaglandidn E2 (PGE 2 ) stimulates osteoblasts to release factors which 
stimulate resorption by osteoclasts. The model equations of Section 3 attempt to incorporate the effects of 
these substances. 


3.0 Model Description 

Current models for bone cell dynamics appearing in the literature of the last decade typically consist 
of three or more equations. The exact number can vary among models. Usually, there is one governing the 
rate of change in the population of the bone forming cells and one governing the rate of change in the 
population of the bone resorbing cells. There may or may not be equations for precursor cells and there 
may be equations for some or none of the biochemicals. In addition, the mathematical descriptions of the 
equations often differ among models. The reason for the difference in the number and form of the 
equations is due in part to the scientific assumptions used to create the equations. Some models are based 
more on biology (Refs. 5 and 7) while some are based more on experimental evidence (Refs. 4 and 21). 

The physical domain in which the cellular dynamics is occurring and upon which the model is based 
varies among models. By physical domain, we mean the type, size, and skeletal location of bone used for 
numerical values of parameters and initial conditions of the model. For example, some models use a 
single BRU (Refs. 4 and 14). Others use an abstract small bone volume (Ref. 5), a section of bone whose 
BRUs are under the same external stimulus (Ref. 7), or a cylindrical section of a long bone (Ref. 9). Only 
a few of the models have an equation governing the rate of change of bone mass in the physical domain. 

In particular, the model of Komorova et al. (Ref. 4) and Moroz et al. (Ref. 14) include a rate equation for 
bone mass as a relative change from the initial value in units of percent. The model of Pivonka et al. 

(Ref. 17) includes an equation for relative change in bone volume in units of percent. Maldonado et al. 
(Ref. 9) use an equation in their model that governs the rate of change in the intracortical radius of a 
cylindrical section of cortical bone to quantity the amount of force induced bone gain. None of the other 
models referenced in the introduction have an equation governing the amount of bone gained or lost or 
the rate of change in the amount of bone. In those models, the gain in the population of active osteoblasts 
(loss in the population of active osteoclasts) or the gain in the population of active osteoclasts (loss in the 
population of active osteoblasts) is related to bone gain or loss, respectively. 

In order to test a model’s ability to simulate bone adaption in the environment of microgravity, the 
model must have a suitable quantification of bone mass at skeletal sites, and it must have a relation 
describing the bone mass as a function of time or an equation describing the rate of change of bone mass 
with time. In addition, the effect of mechanical stimulus and the contribution to bone balance from the 
normal skeletal loads on earth must be included. 

This motivated the consideration of a model that combines the cellular dynamics model equations of 
Lemaire et al. (Ref. 7) with the mechanical stimulus equations in the model of Maldonado et al. (Ref. 9). 
The model in Reference 9 actually uses the cellular dynamics equations in Reference 7 coupled with a 
simple scalar description of a normal earth load on a segment of bone, but not an equation that governs 
bone mass explicitly. Their equation governing intracortical radius appears to have conflicting units 
because the rates of resorption and formation use values from the model of Komorova et al. which are in 
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units that are per cell whereas the rates in the model of Lemaire et al. are in units that are per pM 
concentration. The combined set of equations assemble here avoids conflicting units and attempts to 
model bone matrix lost or gained within a fixed radius. 

According to a study of skeletal data obtained (up to 2007) by the U.S. and Russian space programs, 
trabecular bone loss on a percentage basis is greater than cortical bone loss in certain skeletal sites, but on 
a mass basis, the greatest loss is from cortical bone (Ref. 29). Therefore, the focus herein will be on 
quantifying the change in cortical bone. 

Consider a simple model of a cylindrical section of cortical bone of cortical cross sectional area ,4 and 
length L subject to a mean load F and stress given by 


s = F/A (3.1) 

We will formulate the rate of change of bone area in the cortical region similar to the manner first 
described by Martin (Ref. 11) and further detailed in Reference 12. Figure 2(a) depicts a cross section 
containing secondary osteons or ORUs which are resorbing, refilling, or completed. Pathways of osteons 
may be irregular as well as regular and tunnel almost longitudinally but at a slight angle through the 
cortex (Ref. 12, p. 84). However it will be assumed that osteons and cutting cones of forming osteons are 
perpendicular to the cortical cross section and parallel to the axis of the cylindrical section of cortical 
bone. Therefore, within the cortical cross section there are cross sections of osteons. Figure 2(b) 
represents the cross section of a forming osteon. Martin defines a quantity called intracortical bone 
balance as 


Q - Q b N f - Q c N r 


(3.2) 


where Q B is the mean rate at which bone is replaced during refilling per osteon in mm 2 /day, Q c is the 
mean rate at which bone is resorbed in each osteon during its resorption phase, N F is the number of 
osteons or BRUs forming bone per intracortical cross section A, and N R is the number of BRUs resorbing 
bone per mm 2 cortex. Since resorption begins at the center of a Haversian system and ends at the cement 
line, and formation begins at the cement line and ends at the radius of the canal (Fig. 2(b)), these 
quantities are defined as 


Qb = 




(3.3) 


where R c and R h are the mean radius of an osteon’s cement line and the mean radius of the Haversian 
canal, respectively. Quantities T R and 7> represent the average time it takes for an osteon to complete its 
resorption phase and time to complete refilling phase, respectively. If f a represents the activation 
frequency of a BRU, then in a dynamic situation where the normal remodeling balance is perturbed,^ 
varies with time so that 


N p = \[ +TF f a {s]ds,N R = f 7 * f a (s)ds 

If we assume that the activation frequency/, of a BRU is constant (reasonable for bone in a steady state) 
then 


N F =f a T F , N R =f a T R 


and the times 7> and T R will cancel out of Equation (3.2). 


(3.4) 
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Figure 2. — (a) Sketch of the model. A cylindrical section of cortical bone carries a mean load F. It has a length L and 
cortical area A. Atypical cross section contains Osteon remodeling units which are resorbing (resorption phase), 
refilling (formation phase) or completed, (b) Cross-sectional view of a refilling osteon (formation stage). The radius 
of the cement line is Rq- The radius of the Haversian canal is R^. Formation starts at the cement line, fills inward, 
and stops at the radial distance of the Haversian canal thus creating an area equal to n (R c 2 -Rh 2 )- Prior to 
formation, the resorption phase removed a circular area of bone with radius Rq. 


Actually, under constant^, Equation (3.2) would become a perfectly steady balance (Q = 0) by 
inserting the factor (1 -p) in the second term on the right hand side where 

p = R 2 h /R 2 c (3.5) 

Martin has suggested this alteration of the intracortical balance equation to account for osteons which 
eventually begin to overlap one another, which he has shown to limit porosity in sufficiently remodeled 
bone (Ref. 26). The value of p has been shown to be the upper limit on porosity when osteons are 
randomly located on the cross section (Ref. 26). To further define porosity, we note that the units of Q are 
mm 2 per day or mm 2 /A/day and if A m denotes the area of hard bony matrix, the integral of Q over time is 
the amount of change in A m . Then A-A m , which we denote as A v , is void area consisting of soft tissue. 
Consistent with Martin’s definitions (Refs. 1 1 and 12) the quantity known as bone volume fraction would 
be 


V b = A m /A 


(3.6) 
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and porosity would be 


P v = A v IA (3.7) 

A proposed equation for allowing changes in the normal bone balance Q or dAJdt and a proposed 
equation for the rate of change of calcium in bone fluid are coupled with the cellular dynamics equations 
(Ref. 7) and the mechanical stimulus equations (Ref. 9) to form the model equations. 

4.0 Model Equations 

The tables in Appendix A provide a definition of all the symbols used in this section. A common 
structure can be noted with many of these types of model equations. The rates of change of substances are 
constructed as a combination of accumulation terms and degradation terms. By accumulation 
(degradation) term we mean a term that increases (decreases) the substance rate of change. The system of 
equations here adopts a scheme in which rate coefficients are notated with a plus superscript if they act as 
an accumulation constant and a minus superscript if they act as degradation constant. The subscripts 
associate them with a specific substance. The current model equations are given by the system of 
equations assembled in tabular form below, followed by an explanation for each equation. 


Rate of 
change of 

Equation 

Authorfs) 

Area of hard 
bone tissue 

= QbN f {~) - Q c N r { 1 - p v )(^~) 

at a 0 C 0 

(4.1A) Pennline 

Osteocyte 

concentration 

~~ - r Bo(B - Bq) ~ r Oc(0 ~ Oq) 
at 

(4.2) 

(4.3) Maldonado et al. (Ref. 9) 

(4.4) 

Nitric oxide 
concentration 

AK- no + f - y . r 

^ ' no Jbs no no no 

pge 2 

concentration 

pge _ + r + is r ~ K i T 

£ 'pgeJbs r (no)(pge) n -no ' pge"- pge ' r 1 pge 

OPG 

concentration 

dK 0 R , , . .+ K - K 
U R + l O + r (no)(o)^NO r O K 0 

at tj p 

(4.5) 

(4.6) 

(4.7) Lemaire et al. (Ref. 7) 

(4.8) 

(4.9) 

RANKL 

concentration 

dK L „ . f \ + (k l /k 2 )K 0 + (k i /k 4 )K''^ 

, r L 1 l L n„o)U) K NO r L K L 

dt \ N L E P B J 

Responding 

osteoblast 

concentration 

dB R J p ° 05 d B + „ 

a R tj G £Sr + r ( B )(PSe) K Pge 

Active 

osteoblast 

concentration 

dB = 0 05 d B 
dt e g 

Active osteoclast 
concentration 

^ = d c E R -d A E c C 

Calcium 

concentration 

— = ~ r FB (S/B 0 ) + r BF {Cl C 0 ) + r F + r CF - r FC - r dCa Ca 
at 

(4.10) Pennline 
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The rate of change of A m (t) we propose is given by Equation (4.1 A) where p v is defined in 
Equation (3.7). This accounts for osteon overlap using Martin’s factor discussed in Section 3. Note that 
it’s modeled as a net gain or loss in cortical bone area (defined in Section 3) as a function of the ratio of 
change in cell level to normal cell level for both osteoblasts and osteoclasts. The effects of mechanical 
unloading (and loading) and biochemical imbalances are imposed through changes in cell populations 
caused by these effects. If an equation were used that is similar to the equation governing the rate of 
change of intra-cortical radius in (Ref. 9) it would look like the following 

dA m /dt = Q b N f A - Q c N r {~) + #■+(-£-) - r-(-4a_) (4- IB) 

a 0 L 0 u 0 A m0 

where A m replaces intra-cortical radius, QbN f replaces formation rate in radial direction, QcNr replaces 
resorption rate in radial direction, and and r~ are rate constants. Accumulation due to osteocytes from 
the third term is omitted in Equation (4.1 A) since one could argue that osteocytes are derived from the 
osteoblasts. In comparison to Equations (4.1 A), Equation (4. IB) will be shown to give incorrect results 
for bone loss. 

Since the osteocyte concentration depends on the differentiation of the osteoblast population and the 
osteocytes are subject to cell death or apoptosis, an equation for the rate change of osteocyte 

concentration is given by Equation (4.2), where r RO and r 0c are constant rates of differentiation and cell 
death, respectively. 

Citing references that support their assumptions that osteocytes act as mechanosensors that cause the 
release of Nitrix Oxide (NO) and Prostaglandin E2 (PGE 2 ), Maldonado et al. (Ref. 9) define a nonlinear 
biochemical stimulus function of osteocyte concentration and stress. The biochemical stimulus function is 
given by 


f bs {s,0) 


sO 

(1 + exp(-(r s s + r 0c O))) 


(4.11) 


where r s and r 0c are rate constants. The rates of change of NO and PGE 2 are then defined as a linear 
function of the biochemical stimulus and given by Equations (4.3) and (4.4) where r + 0 , r no , r pge , 

r (no)(pge) > an< i r p g e are rate constants. The quantities I m and I pge represent injection inputs of NO and 
PGE 2 , respectively. 

Lemaire et al. (Ref. 7), describe a chemical kinetics analysis for RANKL, OPG, the binding complex 
of OPG-RANKL, and the RANK-RANKL complex. RANK is held fixed “to reflect the undiminished 
availability of osteoclast precursors”. The four differential equations describing the rate of change of 
these quantities involving the cytokine RANKL and the receptors RANK and OPG are taken from 
(Ref. 6). From those equations Lemaire et al. (Ref. 7) obtain pseudo steady-state expressions for OPG and 
RANKL. In fact, if the steady-state expressions for the complexes OPG-RANKL and RANK-RANKL are 
substituted into the rate equations for RANKL and OPG, the result is Equations (4.5) and (4.6) where 

r (no)(o) is an accumulation constant, r 0 and r [ no )(i) are degradation constants, r L is the rate of R ANKL 

production and elimination, k\lk 2 is the ratio of rates of OPG-RANKL binding/unbinding, and k 2 /k 4 is the 
ratio of the rates of RANK-RANKL binding/unbinding. The constant N L is the maximum number of 
RANKL attached on each cell surface, and r Q is the minimal rate of production of OPG per cell. 
Quantities I 0 and I L are injection rates of OPG and RANKL, respectively. The symbol E P stands for PTH 
receptor occupancy ratio and is a function defined as 
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= Ip(t) + Sp — 

Ip(t) + r p (k 6 /k 5 ) 


(4.12) 


where I P is an injection rate of PTH, S P is the rate of synthesis of systemic PTH, r p is the rate of 

elimination of PTH, and k 6 /k 5 is the ratio of rates of PTH unbinding/binding with its receptor. 

To understand the formation of cell equations, the authors of (Ref. 7) state that their theory is based 
on the assumption that cell “proliferation” (“anti-proliferation”) is proportional (“inversely proportional”) 
to receptor occupancy (Ref. 7, p. 297). The equation governing the rate of change of responding 
osteoblasts is given by Equation (4.7), where d R is a constant representing differentiation rate of 
osteoblast progenitors and E c represents the fraction of occupied TGF-(3 receptors given by 


E G (t) = 


C(t) + 0.05C S 
C(t)+C s 


(4.13) 


The constant Cf represents the value of C to get half differentiation flux. By citing references that 
report that PGE 2 promotes bone formation, Maldonado et al. (Ref. 9) added the last term on the right hand 

side of Equation (4.7). It increases dB R /dt by r (B)( P ge) K pge where f(B)( P ge) is a rate constant. Rates of 

change for the active osteoblasts and active osteoclasts are then described by Equations (4.8) and (4.9) 
where E R = (k^/k A )K L . The constants d B , k B , d c , and d A represent differentiation rate of responding 
osteoblasts, rate of elimination of active osteoblasts, differentiation rate of osteoclast precursors, and rate 
of osteoclast apoptosis caused by TGF-p, respectively. Note that the NO and PGE 2 dynamics is coupled to 
the RANK-RANKL-OPG dynamics to influence the osteoclast and osteoblast interactions through 
Equations (4.5) to (4.7). 

Equation (4.10) is proposed for describing the rate of change of the concentration of calcium in the 
bone fluid and is based on a subset of a model by Jaros et al. (Ref. 3) involving the role of bone in 
calcium homeostasis. The symbols r FB , a- B f, r c f, and r FC represent the rate of flow of calcium from bone 
fluid to bone, rate of flow of calcium from bone to bone fluid, rate of flow of calcium from blood 
capillary to bone fluid and rate of flow of calcium from bone fluid to blood capillary, respectively. The 
constant r F represents the addition of calcium to the bone fluid from bone due to osteocytic osteolysis and 
is affected by the concentration of PTH (Ref. 3). The constant r dCa represents a degradation of calcium 
due to an increased flow of calcium from bone fluid to blood capillary and is set to balance steady state 
for reference values of osteoclasts and osteoblasts. This also accounts for the fact that since an additional 
equation governing the concentration of calcium in the blood plasma is missing; the terms involving r CF 
and r FC are allowed to remain constant. Note that if the full model included an equation governing PTH 
concentration, an additional factor affecting the value of r F would be included. The values for these 
constants are given in mmol per hour for the full skeleton. The units were converted to mmol per day and 
the values were scaled down. See Table 2. 

Equations (4.1) to (4.10) constitute a coupled system of ten 1st order, ordinary differential equations. 
By specifying initial or steady-state values of the variables in Table 1, Equations (4.1) to (4.10) become 
an initial value, first-order system. 


5.0 Simulation and Results 

The quantities Q B and Q c for the formation and resorbtion rates were computed using human rib 
values for the mean Haversian canal and mean cement line radius from Polig and Jee (Ref. 18). An 
average value for activation frequency was computed from a range of rib values referenced by Martin 
(Ref. 11); while an average value of 20 mm 2 for cortical cross-sectional area, A, was chosen, using data 
from Abrams et al. (Ref. 39) which reports cortical area in the cross section of human ribs varying from 
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16.3 to 26.3 mm 2 . The primary intent was to choose numerical values from some published data to be 
near correct orders of magnitude and not to choose cortical values from a specific long bone. They are not 
claimed to be the most accurate currently available. Based on what seemed to be about average normal 
values for the stress on long bones estimated using data from Buckley (Ref. 1, p. 16) on ranges of average 
values of the load in 1 g at various skeletal sites for a 70 kg body weight, a normal load (value of F) on 
earth on the order of 25 N was assumed for the cortical cross sectional area of 20 mm 2 . 

In Equations (4.3) and (4.4), / fo . was replaced by (f bs (s,O)-f hs (25/20,O)) (see Eq. 4.1 1). This was 
necessary in order to be consistent with the definition of the initial cell values for the equations from the 
model by Lemaire et al. (Ref. 7). Their parameter values of reference were selected to represent normal 
bone that is neither gaining nor losing mass, resulting in a system in a state of reference. Therefore, initial 
values from the Lemaire model must be assumed to be from bone with a normal physical load. This 
requires the combined equations to be in steady state when F= 25. Also, in Equations (4.1) to (4.10), we 
replaced the occurrence of each of the variables from Table 1 with the sum of its initial value plus a 
change. For example, B(t) was replaced by B(t)+B 0 , C(t ) was replaced by C(t)+C 0 , and so on. Plots from 
simulations will then show the value of changes from the normal or steady state. 

The method of solution of the equations was as follows. First, following the approach in (Ref. 9), 
Equations (4.3) to (4.6) were solved for their equilibrium expressions by assuming the rates of change are 
zero. This simplification was justified in both (Refs. 7 and 9) by stating that the binding reactions and 
chemical dynamics are much faster than the rest of the dynamics in the model. The expressions for NO, 
PGE 2 , OPG, and RANKL are then substituted in the remaining equations where they appear. This left an 
initial value, first order system of six ODEs, consisting of Equations (4.1), (4.2), (4.7) to (4.10), which 
was numerically solved using a 4th order Runge-Kutta algorithm from Matlab. Unless stated otherwise in 
a specific simulation, the values used for the initial conditions, parameters, and constants are given in 
Tables 2 to 4 and correspond to values in References 7 and 9. 

In reviewing the work of Lemaire et al. (Ref. 7), various bone diseases and possible therapies were 
simulated using their three cell equations. For example, since postmenopausal osteoporosis induced by 
estrogen deficiency is linked to a lack of OPG production (Ref. 27), they ran a simulation by decreasing 
r Q , the secretion rate of OPG by osteoblasts. They then tested a theoretical therapy by reducing k B , the 

rate of elimination of active osteoblasts. In another example, they simulated the catabolic effect of 
continuous injection of PTH by running a simulation with I P set to 1000 pM/day. 

To check the agreement of the computations from the model Equations (4.1) to (4.10) with those of 
Lemaire et al. and the ability of the model to simulate bone disease on earth with no change in the normal 
force, a simulation with I P set to 1000 pM/day was performed. Figure 3 shows the results wherein the 
normal load of (F(t) = 25) was used. This reduced Equations (4.7) to (4.9) to the same three cell equations 
employed in (Ref. 7) to give results for the osteoclasts and osteoblasts curves identical to those displayed 
in (Ref. 7, p. 300). The difference here is that simulation results showing the change in the area of hard 
bony tissue in the cortical cross section (lower right; Fig. 3) and the change in bone fluid calcium (lower 
left; Fig. 3) are displayed, whereas they are not displayed in (Ref. 7). 

Note that the curve resulting from the use of (4. IB) shows an increase in bone area during the 
continuous injection period which is inconsistent with experimental and clinical observations (Ref. 30), 
whereas the response resulting from Equation (4.1 A) shows the expected decrease in bone area. Actual 
values of the result from Equation (4. IB) were reduced by a factor of 10 _1 to enable comparing the 
qualitative difference on the same plot. 

In reviewing the use of their force-induced bone growth model, Maldonado et al. (Ref. 9) compared 
the difference between the effects of bone disease on the intracortical radius with physical activity and the 
effects without physical activity. They were able to show that physical activity is beneficial. For example, 
in the first test of their model they simulate the effects of an increased load on the intracortical radius. 
They ran a 350 day simulation during which the force F was increased on day 50 to mimic physical 
activity which increased their assumed average value of F to a value higher than the assumed normal 
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Figure 3. — Plots showing changes in osteoclasts, osteoblasts, bone fluid calcium, and bone area from 
a continuous injection (1000 pM/day) of PTH started on day 20 and then stopped on day 80. 


value over the next 100 days. It was further increased on day 1 50. On day 250 the force was dropped to a 
value lower than the assumed normal value, to mimic a decrease in physical activity over next 100 days. 
The results show an increase in the radius between days 50 and 150 and another increase between days 
150 and 250. The last one hundred days showed a decrease in the radius to a level slightly below the 
equilibrium value. A figure of those results can be seen in Reference 9. 

For purpose of this study, a 350 day simulation of the model using Equation (4.1 A) was conducted, 
where the load (value of F) was first decreased to zero on day 50. This would correspond to unloading the 
cylindrical section of bone upon entrance into the environment of microgravity on day 50. Figure 4 shows 
the results of such a simulation which also includes the bone area plot obtained using Equation (4. IB). 

Once again, the response (change in bone area) resulting from the use of Equation (4.1 A) predicts 
expected qualitative results in bone loss and gain. There is a continued loss of bone when the load is 
decreased to zero on day 50. On day 150 the average load is increased over the next 100 days, as the 
cumulative effects of physical activity is accounted for by increasing F to 20, and there is a decrease in 
the rate of loss of bone area, as shown. Corresponding to a return to 1 g on day 250, the load is increased 
above the normal load and bone loss begins to reverse. 

On first glance, the response resulting from the use of Equation (4. IB) also indicates expected 
qualitative results in bone loss and gain in the following sense. Bone area in the cross section is lost in the 
interval of time beginning with a zero load. Raising the average load to 20 over an interval of time 
beginning on day 150 recovers some of the loss. Leaving the environment of a reduced load at day 250 
and increasing the average load to 30 with physical activity recovers the lost bone area and then some. 

The shape of these curves is similar to the shape of the curves obtained in Reference 9 for bone growth 
induced by increasing the average force. However, it is not a good model of bone loss in the environment 
of microgravity where bone loss can continue over long periods of time, since the bone curve appears to 
level off or almost approach an asymptotic value in the unloaded interval (F = 0) without any 
countermeasure. 
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Figure 4. — Plots showing the changes in osteoclasts, osteoblasts, bone fluid calcium, and bone area 
obtained by setting F - 0 on day 50, increasing Fto 20 on day 150, and then increasing the average 
load F to 30 on day 250. 


Thus, the results of the present simulation demonstrate the following. The model Equation (4.1 A) for 
rate of change of bone, which is more biologically accurate than Equation (4. IB), is able to more 
correctly simulate expected bone loss in microgravity which, in this model, means by reducing the normal 
value of F {F = 25) assumed on earth. 

Biochemical imbalances under the conditions of microgravity may similarly cause disruptions in the 
bone remodeling balance in addition to those disruptions explained or modeled strictly by the mechanics 
of skeletal unloading. Lemaire et al. (Ref. 7) remark that osteoblast progenitors express TGF-P receptors 
so that TGF-P, produced by osteoclasts during bone resorption, leads to differentiation of the progenitors 
into responding osteoblasts. Using that knowledge, they associate bone loss, due to senescence on earth, 
with a decrease in the production of TGF-p. However, it is not clear what numerical values they modify 
to simulate bone loss due to senescence. We associated the consequence of senescence with one of the 
possible causes of old age osteoporosis, a decrease in the rate of osteoclast apoptosis due to TGF-p. 

Figure 5 shows the result of such a simulation of the model equations (without any skeletal unloading) 
using Equation (4.1 A) for bone balance. Notice that when the osteoclast death rate due to TGF-P is 
reduced from 0.7 to 0.5 in the interval between day 20 and day 80, the loss of bone is similar to the loss in 
Figure 4 during the interval between day 50 and day 150. Bone area recovers and returns to the balanced 
state when the osteoclast death rate due to TGF-P is returned to its normal value. 

It is of interest to examine the results of a simulation of this model that combines unloading with a 
reduction in the rate of osteoclast apoptosis due to TGF-P as a possible combination of bone remodeling 
imbalances in the environment of microgravity. Figure 6 shows the results of such a simulation using 
Equation (4.1 A) for bone balance. 
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Figure 5. — Plots of bone fluid calcium in the section and bone area in the cross section. The value of 
d/\, the osteoclast death rate due to TGF-p, was reduced from 0.7 to 0.5 on day 20 and increased 
back to 0.7 on day 80. The simulation used Equation (4.1A). 




Figure 6. — Plots of changes in bone fluid calcium and bone area. Between day 50 and day 250 d/\, was 
reduced from 0.7 to 0.5. The load was set to 0 on day 50 and to 50.0 on day 150. On day 250, the 
load remained at 50.0 and dfy was restored to 0.7. The bone balance Equation (4.1A) was used. 


Observe that the combined effects of reducing the rate of osteoclast apoptosis and unloading make it 
harder to recover bone loss with say physical activity that increases the average load. Increasing the 
average value of F above the normal load to 50.0 over the interval between day 150 and 250 still only 
showed a decrease in the rate of bone loss whereas without the added effect, bone loss would have been 
reversed. Other therapies and countermeasures besides increasing F could also be tested. However, we 
defer additional studies until a better model is developed, especially a model with a more robust 
description of the mechanical load on bone and mechanotransduction. 


Conclusions 

The simulations carried out with the simple model in this study show that mathematical formulations 
of the bone remodeling function could eventually be used for predicating bone loss in space. It is 
important that the model equations have a steady state or semi-steady state accurately perturbed by well- 
defined metabolic and biomechanical scenarios. However, to have a model with sufficient high fidelity to 
be used for application in the development of DA, or any other realistic simulations, much more research 
and development is needed. How bone loss or gain is quantified needs to be more accurately linked to the 
biological structure and remodeling phase of bone. Theories on the exact biomechanical stimulus for 
strain adaptation and mechanotransduction need to be more completely linked to the metabolic process. 
To be more specific, consider some of the limitations and flaws of the model in this study. 

Although qualitative results look interesting, the modeling is still too abstract to be used for 
quantitative purposes, and this is true for other models in the literature. There is not enough dimensional 
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detail to be able to associate a collection of remodeling units under the same external stimulus with a 
specific location and size of bone segment. Our model does not include the complete set of cell types and 
biochemicals. For example, the interactions of phosphate ions are missing. The remodeling of trabecular 
bone is not addressed. The model here (and most other models) alters the resorption or formation rate 
based only on a change in the concentration of cells, and does not include a relation to changes in 
resorption depth or direction that Parfitt (Ref. 15, p. 73) associates with certain bone loss conditions. 

More research in needed in order to fully understand how the applied loads induce osteogenic response by 
various cell types in different regions (Ref. 16). Work is needed in designing experiments to quantify all 
model parameters and to insure unit and dimensional consistency, developing a benchmark for meeting 
requirements for realistic simulations, and developing validation tests. 
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Appendix A. — Tables 


TABLE 1.— NOMENCLATURE 


Symbol 

Definition 

Unit 

4ni 0 

Area of bone 

mm 2 

B(t) 

Concentration of osteoblasts 

pM 

B R {t ) 

Concentration of responding osteoblasts 

pM 

dt) 

Concentration of osteoclasts 

pM 

0(0 

Concentration of osteocytes 

pM 

Ca(t) 

Calcium concentration 

mmol 

K,M 

Concentration of RANKL 

pM 

K„ 0 (t) 

Concentration of nitric oxide 

pM 

K a (t ) 

Concentration of osteoprotegerin (OPG) 

pM 

K p g e(t) 

Concentration of prostaglandin E 2 (PGE 2 ) 

pM 

hiO 

Injection rate of RANKL 

pM/day 

ho(t) 

Injection rate of nitric oxide 

pM/day 

hit) 

Injection rate of OPG 

pM/day 

hit) 

Injection rate of PTH 

pM/day 

hgeit) 

Injection rate of PGE 2 

pM day' 1 


TABLE 2.— INITIAL VALUES 


Symbol 

Value 

Definition 

A m0 

18 mm 2 

Initial value bone area 

Bo 

0.0007282 pM 

Initial value active osteoblasts 

Bro 

0.0007734 pM 

Initial value responding osteoblasts 

Co 

0.0009127 pM 

Initial value active osteoclasts 

C 

0.0073 pM 

Initial value osteocytes 

Ca 0 

0.696 mmol 

Initial value calcium (Bone fluid) 
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TABLE 3.— PARAMETER VALUES 



Value 

Definition 


0. 173x 10 -3 mm 2 day -1 

Rate of formation of bone area by osteoblasts 

QcNr 


Rate of resorption of bone area by osteoclasts 

r m 

1.0 mm 2 day” 1 

Rate of increase of bone area due to osteocytes 

r m 

1.0 day -1 

Rate of degradation of bone area 

r BO 

0.1 day -1 

Osteocyte production rate from differentiation of osteoblasts 

r Oc 

1 day -1 

Osteocyte degradation rate 

r s 

1 mm 2 /N 

Stress influence rate 

r Oc 

1.0 pM" 1 

Osteocyte influence rate 

r no 

2.0x 10 4 mtn 2 day/N 

Rate of release of NO 

r no 

l.OxlO 3 day -1 

NO elimination rate 

r + 

Pge 

] 00 mm 2 day/N 

Rate of release of PGE 2 due to biochemical stimulus 

r + 

( «o)(pge ) 


PGE 2 rate increased by NO 

r~ 

Pge 

l.OxlO 2 day 1 

PGE 2 elimination rate 

r _0 

2xl0 5 pmol day '/pmol cells 

Minimal rate of production of OPG per cell 

r (no)io) 

10 day 1 

OPG rate increase by NO 

r o 

0.35 day 1 

Rate of elimination of OPG 

r L 

10 3 pM day ' 

Rate of RANKL production and elimination 

r (no)(l) 

1 00 day 1 

RANKL rate decreased by NO 

r p 

86 day”' 

Rate of elimination of PTH 

r (B)(pge) 

1.0x10 4 day ’ 

B r rate increased by PGE 2 
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TABLE 4.— PARAMETER VALUES 


B3® 

Value 

Definition 

c s 

5 X 10 -3 pM 

Value of Cto get half differentiation flux 

d A 

0.7 day -1 

Rate of osteoclast apoptosis caused by TGF p 

d B 

0.7 day -1 

Differentiation rate of responding osteoblasts 

d c 

IBS; 

Differentiation rate of osteoclast precursors 

d R 


Differentiation rate of osteoblast progenitors 

K 

10 pM 

Fixed concentration of RANK 

k\ 

gftjfjpfi .v 

Rate of OPG-RANKL unbinding 

^2 

10 day 1 

Rate of OPG-RANKL unbinding 

k 3 


Rate of RANK-RANKL binding 

k 4 

1.7xl0“ 2 day' 1 

Rate of RANK-RANKL unbinding 

k 5 

ggjggpggl^gBI 

Rate of PTFI binding with its receptor 


3 day 1 

Rate of PTFI unbinding 

n l 

3><10 6 pmol/pmol cells 

Maximum number of RANKL attached on each cell surface 

k B 

0.1 89 day" 1 

Rate of elimination of active osteoblasts 

r FB 


Rate of flow of calcium from bone fluid to bone 

r BF 


Rate of flow of calcium from bone to bone fluid 

r CF 


Rate of flow of calcium from blood capillary to bone fluid 

r FC 


Rate of flow of calcium from bone fluid to blood capillary 

r F 


Rate of accumulation of calcium in bone fluid due to osteocytic osteolysis 

r dCa 

0.1 day 1 

Calcium degradation rate 
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Appendix B. — Substance Definitions 


Unless specifically referenced, definitions were taken from Wikipedia, The Free Encyclopedia, 

http://en.wikipedia.org/wiki/. 

1 . Receptor. — A protein molecule, embedded in either the plasma membrane or cytoplasm of a cell, to 
which a mobile signaling molecule may attach. 

2. Ligand. — A molecule that binds to a receptor. 

3. Receptor Activator of Nuclear Factor k B (RANK). — Membrane protein which is expressed on the 
surface of osteoclasts and is involved in the activation of osteoclasts upon ligand binding. 

4. Receptor Activator for Nuclear Factor k B Ligand (RANKL). — A surface bound molecule expressed 
on osteoblasts that stimulates osteoclast differentiation from preosteoclasts through binding with 
RANK. 

5. Parathyroid Hormone (PTH). — A hormone secreted by parathyroid gland that enhances release of 
calcium from bone and indirectly stimulates osteoclasts by binding to osteoblasts which in turn 
stimulates osteoblasts to increase expression of RANKL. 

6. Osteoprotegerin (OPG). — A decoy receptor for RANKL that prevents osteoclast maturation by 
blocking the RANKL-RANK interaction. 

7. Transforming Growth Factor beta (TGF-P). — Controls cellular differentiation and other functions in 
most cells and induces apoptosis in numerous cell types. TGF-P is released by osteoclasts during bone 
resorption. It has the ability stimulate addition in and activity in the preosteoblast population but it 
can also inhibit final maturation into active osteoblasts (Refs. 35 and 36). It’s also known to induce 
osteoclast apoptosis (Ref. 37). 

8. Nitric Oxide (NO). — A gaseous signaling molecule produced by osteoblasts which directly regulates 
osteoclastic activity (Ref. 34). 

9. Prostaglandin E2 (PGE 2 ). — A substance that stimulates osteoblasts to release factors which stimulate 
bone resorption by osteoclasts. 
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